Physics of the interior of a spherical, charged black hole with a scalar field 
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We analyse the physics of nonlinear gravitational processes inside a spherical charged black hole 
perturbed by a self-gravitating massless scalar field. For this purpose we created an appropriate 
numerical code. Throughout the paper, in addition to investigation of the properties of the mathe- 
matical singularities where some curvature scalars are equal to infinity, we analyse the properties of 
the physical singularities where the Kretschmann curvature scalar is equal to the planckian value. 
Using a homogeneous approximation we analyse the properties of the spacetime near a spacelike 
singularity in spacetimes influenced by different matter contents namely a scalar fleld, pressureless 
dust and matter with ultrarelativistic isotropic pressure. We also carry out full nonlinear analyses 
of the scalar field and geometry of spacetime inside black holes by means of an appropriate numer- 
ical code with adaptive mesh refinement capabilities. We use this code to investigate the nonlinear 
effects of gravitational focusing, mass inflation, matter squeeze, and these effects dependence on the 
initial boundary conditions. It is demonstrated that the position of the physical singularity inside a 
black hole is quite different from the positions of the mathematical singularities. In the case of the 
existence of a strong outgoing flux of the scalar field inside a black hole it is possible to have the 
existence of two null singularities and one central r = singularity simultaneously. 

PACS numbers: 04.70.Bw, 04.20.Dw 
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I. INTRODUCTION 

The problems of the internal structure of black holes 
are a real great challenge and has been the subject of 
very active analytical and numerical researches durinj 
the last decades [u p, II , 1 1. [a|a , H , [a 

been a great progress in these researches in the last iew 
years and we now know many important properties of 
the realistic black hole's interior, but some details and 
crucial problems are still the subject of much debate. 

Many important results have been obtained under sim- 
plifying assumptions. One of the most widely used test 
toy models is a spherical, charged, non-rotating black 
hole, nonlinearly perturbed by a minimally coupled and 
self-gravitating massless, uncharged, scalar field. While 
this toy model is not very realistic, it share many prop- 
erties, including causal structure, with the more realistic 
rotating black holes (e.g. page 5 and which is 
why it is believed that insights into this model may give 
us important understandings about rotating black holes. 

The purpose of this paper is to continue the analysis 
of the physical processes in the interior of black holes in 
the framework of this toy model. 

Inside a black hole the main sights are the singularity. 
A number of rigorous theorems (see references in |2g) 
imply that singularities in the structure of spacetime de- 



velops inside black holes. Unfortunately these theorems 
tell us practically nothing about the locations and the 
nature of the singularities. It has been found that in 
principle two types of singularities can exist inside black 
holes corresponding to the toy model: A strong spacelike 
singularity and a weak null singularity (instead of the in- 
ner horizon of a Reissner-Nordstrom black hole). Prob- 
ably both types of singularities can exist simultaneously 
in the same black hole and probably it is possible to have 
cases where o nly the strong spacelike singularity exists. 
In the works |23 and references therein, some physical 
and geometrical properties of the singularities have been 
investigated. Numerical simulations of the fully nonlin- 
ear evolution of the scalar field and the geometry inside 
the spherical charged black hole has been carried out in 
(Id .14. .15 . .2Sj . 

Near the strong spacelike singularity, one can use a 
homogeneous approximation in which it is supposed that 
temporal gradients are much greater then the spatial gra- 
dients. Hence it may be assumed that all processes near 
the singularity depend on the time coordinate only. This 
uniform approximation has been used in (page 212) 
and |llL I23 to gain important new knowledge about the 
processes near the singularity. We will use the same ho- 
mogeneous approximation to extend these analyses and 
clarify some fundamental physical processes near space- 
like singularities under the influence of three different 
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matter contents, namely for the case of pressureless dust, 
a massless scalar field and matter with ultrarelativistic 
isotropic pressure. This investigation is done by means 
of a suitable numerical code which we develop for this 
purpose. 

Subsequently we will study the nonlinear processes in 
large regions inside the toy model black hole, not just 
limited to the homogeneous approximation near the sin- 
gularity. This will be done by using a stable and sec- 
ond order accurate, numerical code with adaptive mesh 
refinement capabilities. We will perturb the black hole 
with initial infalling scalar fields of different forms and 
strengths to further investigate the behaviour of the sin- 
gularities and physical processes near them. We will also 
investigate the influence of outgoing scalar fluxes on the 
interior regions and the singularity. Such outgoing fluxes 
will unavoidably appear as a result of the scattering of 
ingoing scalar fleld flux by the curvature of the space- 
time and will also be emitted from the surface of a star 
collapsing to a black hole. 

Today it is widely believed that in the singularity of a 
realistic black hole, the curvature of the spacetime tends 
to infinity. Close to the singularity, where the curvature 
approaches the Planck value ((^) ~ 1-5 • 10^^^ cm"'' 
pjQj). classical General Relativity is not applicable. There 
is not yet a final version of the quantum theory of gravity, 
thus any extension of the discussion of physics in this 
region would be highly speculative and we will consider 
these regions as singularities from the classical point of 
view throughout the paper. 

The paper is organised as follows; In section ^ 
we present our model of the spherically symmetric, 
charged black hole. In section IIIII we dicuss the 
mass function and some important nonlinear effects 
which are fundamental for understanding the physical 
processes inside black holes. In section IIVI we use a 
homogeneous approximation to analyse the (spacelike) 
singularity for three different matter contents: dust 
with zero pressure, a massless scalar field and matter 
with relativistic isotropic pressure. In section we 
analyse the full nonlinear equations of the model from 
section ^ using a numerical approach. Finally we 
summarize our conclusions in section I VII Details of 
the numerical code used to obtain the results in sec- 
tionlVland analysis of it are given in appendices IXI and IbI 



II. THE MODEL 

We wish to study the geometry inside a spherically 
symmetric black hole with a fixed electrical charge q (i.e. 
Reissner-Nordstrom metric), which is nonlinearly per- 
turbed by a selfgravitating, minimally coupled, massless 
scalar field. While astrophysical black holes are more 
likely to be described by the Kerr metric, it is believed 
that this toy model captures the essential physics, 
since the causal and horizon structures of the Reissner- 



Nordstrom and Kerr black holes are known to be very 
similar |2[ (page 5). However, it is much simpler to make 
a numerical model of the toy model since this can be 
simulated in a two-dimensional spacetime. In construct- 
ing the toy model, we follow here the approach of Burko 
and Ori p!4 . ll5ll28| who have done similar investigations. 



A. Field equations 

In spherical symmetry, the general line element in dou- 
ble null-coordinates can be written as: 



(1) 



where dVl^ = dO^ + s\v?{9)d(jP' is the line element on the 
unit two-sphere and r is a function of the null coordinates 
u and V (in- and outgoing respectively). 

With this metric the non-zero components of the Ein- 
stein tensor are: 
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= -2e ^"rir^uv +r(7^uv) 
= -2 e-^" r sin2(6') (r,„„ +ra,uv) (2e) 

The energy-momentum tensor can be written as a sum 
of contributions from electromagnetic and scalar fields: 



rp rpS I rpe7n 



(3) 



The energy-momentum tensor of a massless scalar field 
$ is 113: 

T^. = ^ - ^5M^5"'''i>,a$,/3) (4) 

whose non-zero components for the metric Q are: 
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sm^{e)e-'"'<S>^u'i>,v (5d) 



The energy-momentum tensor of an electric field in 
spherical symmetry and null coordinates is |3d] : 



rpem p pa p pfii. 



(6) 
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whose non-zero components for the metric Q are: 
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From the Einstein and energy-momentum tensors we 
can write up the Einstein equations, G^i, = 8 tt T^^ (with 
c = 1, G = 1), governing the spacetime. The u — u, v — v, 
u — V and 9 — 9 components of the Einstein equations 
respectively are: 



J',™ - 2r,„CT„ -l-r ($,„) =0 
r_„^ - 2 r^i,a „ + r ($,1,)^ = 



r 2r 



2a / 2 
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(8) 
(9) 

(10) 
+ - 0(11) 



Lastly, the scalar field must satisfy the Gordon-Klein 
equation (note that Gordon-Klein is a consequence of the 
Einstein equations for the scalar field (30,]), V^V^$ — 0, 
which in the metric becomes: 



r 



(12) 



Equations (|10ll - (|12|l are evolution equations which 
are supplemented by the two constraint equations 
and lO . It is noted that none of these equations depends 
on the scalar field $ itself, but only on the derivatives 
of <&, i.e. the derivative of the scalar field is a physical 
quantity, while the absolute value of the scalar field itself 
is not. Specifically we note the Tuu — {^.u)'^/4:Tt and 
Tvv — (*i',D)^/47r components of the energy-momentum 
tensor which are part of the constraint equations. 
Physically T^u and Tyy represents the flux of the 
scalar field through a surface of constant v and u 
respectively. These fiuxes will play an important role 
in our interpretation of the numerical results in sectionlVl 



B. Initial value problem 

We wish to numerically evolve the unknown functions 
r{u,v),a{u,v) and $(?/,«) throughout some computa- 
tional domain. We do this by following the approach of 
P,0, 15, 28] to numerically integrate the three evolution 
equations ifTUI) - l(T^ . These equations form a well-posed 
initial value problem in which we can specify initial values 
of the unknowns on two initial null segments, namely an 
outgoing {u = uo = constant) and an ingoing (v = vq = 
constant) segment. We impose the constraint equation 



^ and Q on the initial segments. Consistency of the 
evolving fields with the constraint equations is then en- 
sured via the contracted Bianchi identities but we 
use the constraint equations throughout the domain of 
integration to check the accuracy of the numerical simu- 
lation. 

On the initial null segments, the constraint equations 
reduces the number of unknowns by one on v = vq and 
u = uq respectively. The remaining two unknowns ex- 
presses only one degree of physical freedom, while the 
other unknown expresses the gauge freedom associated 
with the transformation u —^ u{u),v — *■ v{v), (the line 
element and eqs. (|5|)- (I12|I are invariant to such a 
transformation) . We choose a standard gauge in which r 
is linear in v and u on the initial null segments. Specifi- 
cally we choose: 



r{uo,v) 



r{u,vo) = ro -|- tir„o 



(13) 



We also choose that the outgoing segment should run 
along Uo ~ 0. 

We can now use © and © to find: 

a^y{uo,v) = (r^yy +r($,„)^^ = ^ ($,„)^ (14a) 



2r. 



" ($.,)^ 



(14b) 

which can can be readily integrated to find cr(u, v) on the 
initial null segments if $ and the constants rg , r„o and 
cr(uo,fo) are specified on these. 

Following [ll m m we choose a{ua,vo) = 
and ro = 5. The parameter r„o, can be related to the 
initial mass and charge of the black hole via the mass 
function (the mass function is further discussed in section 
llll|) which in the metric ^ has the form: 



m{u, v) 



1 



2e2 



(15) 



which in our choice of gauge, at the point of intersection 
of the initial null segments, take the form: 



mo = m(uo,wo) 



^•0 



1 



4r„ 



hence r„o can be determined by ro, too and q as: 

^2 



1 /2 / q 



- 1 



(16) 



(17) 



Hence, by specifying a distribution of the scalar field 
$ on the initial null segments, choosing a gauge and ini- 
tial charge and mass of the black hole we can specifiy 
complete initial conditions on the initial null segnemts. 
Using a numerical code (described in appendix ^ we 
can then use the evolution equations, eqs. H10|) - (|12|) to 
evolve the unknown functions throughout the computa- 
tional domain. 
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III. NONLINEAR EFFECTS; INTERNAL MASS 
FUNCTION 

We will in the next sections consider the evolution of 
the scalar field together with the geometry of the inte- 
rior of a black hole. This evolution is highly nonlinear. 
One of the main parameters of this evolution is the mass 
function which represents the total effective mass in a 
sphere of radius r{u, v) 0,11^1^. We give here different 
expressions for the mass function, which emphasize its 
different characteristics. 

In the metric: 

ds^ — gttdt^ + grrdr^ + r^dVl^ 
dfl^ = d0^ + sin^ ed(t>^ (18) 

the mass function can be written in the following forms: 

"1 = ^ (l + ^ - 9rr^ (19) 

or 

m = 47r / T*r^dr + mo (20) 

Jri 

In the metric 

ds^ = -a^dudv + r^dn"^ (21) 
it has the form [Tll28t: 

or (for the scalar field, fl3 |: 

m^™ = 2-$2^$2^ _ , (^1 _ _ + ^ j (23) 

There are two important physical processes which can 
lead to a nonlinear change of the mass parameter: 

1. The mass m inside a sphere can change because of 
the work of pressure forces on the surface of the 
sphere. A clear manifestation of this squeeze effect 
is the change of the mass of a spherical volume in 
a homogeneous model of the Universe filled with 
relativistic gas (see [13 1 page 13). In section Hvl 
we will consider another example, namely for the 
case of the imitation of the interior of a black hole. 
For the description of this process, it is most ap- 
propriate to use the form H2U|) for the mass func- 
tion. Remember that inside the event horizon, r is 
a time-like coordinate. 

2. Mass inflation This process inside the black 
hole exists if near the Cauchy horizon (CH) there 
are simultaneous ingoing and outgoing fluxes of a 
massless field (for example scalar field). Actually 



the existence of the outgoing flux together with the 
ingoing is inevitable because of backscattering of 
part of the ingoing flux by the spacetime curvature. 
The simplest exact model of the mass inflation pro- 
cess inside of a black hole has been constructed by 
Ori 0. For the description of this process it is 
most useful to use formula H23|l from which it can 
be seen that evolution of m with both u and v is 
possible only if there are both ^\ and ^\ fluxes 
simultaneously. 

One can often observe the simultaneous manifestation of 
both these processes. 

Another important nonlinear effect is the focusing ef- 
fect by the gravity of beams of opposite fluxes of radi- 
ation. A particular manifestation of the focusing effect 
is the contraction of the CH under the gravity of trans- 
verse irradiation by the outgoing radiation. Eventually 
the CH singularity shrinks down to a point-like size and 
meets a central (probably spacelike) singularity r = 0. It 
should be mentioned that it is incorrect to say that the 
CH singularity is transformed into a r = spacelike sin- 
gularity, because the formation of the r = singularity 
is causally absolutely independent from the formation of 
the CH singularity. 

We want to mention that it is possible, in principle, to 
have the situation when the mass function depends on 
only one null coordinate, say w, while it does not depend 
on the other u coordinate. This situation is described by 
a charged Vaidya solution [s^. In this solution there is 
an effect of a linear change of m, because of an ingoing 
lightlike radial flux of energy into the black hole (with- 
out any scattering of this radiation by a curvature of the 
spacetime). Of course this effect is compatible with for- 
mula 

IV. HOMOGENEOUS APPROXIMATION 

In the close vicinity of the spacelike singularity of a 
black hole all processes, as a rule, have high temporal 
gradients, much higher than the spatial gradients along 
the singularity, and the processes depend on a very re- 
stricted space region. So for clarification of some physical 
processes one can use a homogeneous approximation and 
assume that all processes depend on the time coordinate 
only. This approach has been proposed by Burko 0, 
and we will use it at the beginning of our analysis to clar- 
ify some main properties of the singularity before coming 
to the full analysis of the spherical model in section Ivl 

1. Leading order analysis 

The general homogeneous, spherically symmetric line 
element has the form: 

ds^ = gtt{r)dt'^ + grr{r)dr^ + r^dQ'^ 
dfl^ = d0^ + sin^ 0d(l)^ (24) 
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Inside a black hole in the region between the event hori- 
zon and the Cauchy horizon (or the spacelike singularity) 
r is timelike and t is spacelike. To describe the contrac- 
tion of the CH, we should thus consider the variation of 
the time coordinate r from bigger to smaller values. The 
r — r, t — t and — components of the Einstein equations 
(with c = 1, G = 1) are then given by: 



9tt - 9rr 9tt + rg'u 



"T grr 9tt 



2 / 

r2 g^"^ 



= 87r(T;-H^;) 
= Svr (T* + i?*) 



(25) 
(26) 



Argrr^ gtt 



h-^{9tt[2grr{gtt + rg't't) - {rg'rrg'u)] 



-2gt^ g'„ - rgrr g'u } = MT'e + E^) (27) 

where the primes denotes differentiation with respect to 
r (the (j) ~ 4' component of the Einstein equations again 
yields equation ^^). The tensor E represents here con- 
tribution from a free electric field corresponding to a 
charge g, which we will assume to be constant: 



e: = El 



-E' 



(28) 



while the tensor T represent contributions from other 
matter contents. 

To clarify the meaning of different processes we will 
consider three different physical matter contents (in ad- 
dition to the electric field) with different equations of 
state. Namely, we will consider: 

A) Dust 

B) A massless scalar field 

C) Ultrarelativistic gas 

From the Einstein equations one can find the following 
expressions for the non-zero components of T for these 
matter contents: 



A) Dust (vifith pressure P = 0) 

/ gtt^init 



-e = eo 



V 9tt 

where eo,gtt,imt and ri„it are constants 
B) Massless scalar field j^: 
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StT gu,init'''init 



(30a) 
(30b) 

(30c) 

2 

(30d) 
(30e) 



where d is a constant which were used in . 

C) Ultrarelativistic gas (with isotropic pressure 
P = |, e being matter density): 

(31a) 
(31b) 

(31c) 



t: = -e 



3 

9tt,init 

gtt 



^ '^init ^ 



Substitution of l(^ - (p?T|l into and enables us 
to find the unknown functions grr — grr(r) and gtt = 
gtt (r) and hence solve the problem for each of the three 
cases. Equation (|?7jl can be used as a control of the 
calculations. 

Formally metric H24|) corresponds to a metric of a spe- 
cial class of the "homogeneous cosmological models" con- 
sidered by Zeldovich and Novikov [331 (page 535), Gr- 
ishchuk f34!] and others. 



A. Dust, P = 

Let us start the discussion with the simplest case, 
namely the case of dust with pressure P = 0. 

In this case it is possible to have a singularity at r = 
Tsing 7^ with rcH < r sing < TEH, whcrc rcH and teh 
are the positions of the Cauchy Horizon and the Event 
Horizon in the absence of dust. 

Let us consider the leading order terms in a series ex- 
pansion for the metric functions and leading order terms 
in the Einstein equations, near the singularity. Close to 
the singularity, where gtt — > 0, a leading order expansion 
of eqs. + gives us: 



dgt 
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{grr)s 



£0 



dx gttTs 



(32) 



where {grr)s 



[rsing) and gtt. 



gtt,init ij'init ) 



and where we assume that gtt — Ax"^ and x = r 



Also a are constants and {grr)s 
value of grr at the singularity r = 
From H32I) one find: 



a 



grr (rsing) is the 



(33) 



which leads in turn to (remember that {grr)sing is nega- 
tive for rcH < r < veh)- 

_ 47r {grr)sing Ep rf^it ^gtt,init ^^^^ 



' Sing 



A 



Using the proper time t : dr — \/\grr\dr we have for the 
vicinity of the singularity: 

2 



gtt oc r 
r = const., 
T = at the singularity 



(35) 
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FIG. 1: Metric coefficients (plot a) and Kretschmann scalar (plot b) versus r for the case of dust with eo= 0.03 



For the Kretschmann scalar K = i?**^''™ we have 



K 



12 



(36) 



Thus this spacelike singularity does not correspond to 
r = 



1. Numerical analysis 

To understand the behaviours of the model H24|l as 
functions of the parameters of the model we perform use 
a simple numerical code to numerically solve H25() + l|26() 
substituting (|29|l for the stress-energy tensor. According 
to the remark in the introduction we consider the region 
with K ~ Kpianck as a physical singularity and will con- 
sider only the region with r > rc, where Tc corresponds to 
the critical value of r at which the Kretschmann scalar is 
equal to the planckian value. We will take this restriction 
into account in all our subsequent analyses. 

We note that for the case of dust P — there are no 
nonlinear effects causing an increase of the mass function 
771. This is seen from eq. H2Uf) because — P = 0. 

To analyse the change of Tc with variation of the matter 
contents we numericallv integrate eas. (|25|l . H2()l) . (|29|l . As 
initial values we use rinit = 0.95 • veh ~ 1-25 and set 
gtt.init and grr,init equal to their values at rinit for the 
Reissner-Nordstrom solution (with initial mass rriQ — 1 
and charge q = 0.95) and vary the initial matter density 
eo- 



Figure 1(a) shows an example of the variation of the 
metric functions with r for the case cq = 0.03. It is 
seen that gtt ^ as r ^ Tc w 1.149. As gu{r) 0, 
the density and curvature increases rapidly, which can 



easily be un derstood from eq. H29|) . This is indicated 
in fig. l(b)| which shows the variation of K{r) for the 
same case. The line in this figure does not visibly reach 
K — Kpianck ~ 1.5 • 10^^^, however this is solely due to 
limitations in numerical resolution because of the catas- 
trophic blowup of K{r) as indicated by the vertical line 
in the figure. 

The dependence of on the initial matter density eo 
can be seen in fig. |2] Near the mathematical singularity, 
fsing, the scalar K increases very rapidly with decreas- 
ing r, so approximately r^ing ~ rc (physical singularity). 
As one can see from the figure, for the case of dust, rc 
decreases with decreasing eo until rc — > rcH at eo ^ 0. 
This behaviour is easily understood: for smaller matter 
contents it takes a longer time to compress the dust to the 
critical density at rc- On the other hand, in the Reissner- 
Nordstrom solution without additional matter, the vol- 
ume of the uniform reference frame H24|l tends to zero 
when r rcH (because gu —^0). So when eo ~> and 
the behaviours of the solutions are close to the Reissner- 
Nordstrom solution, the matter density of dust must tend 
to infinity when the volume tends to zero at r — > rcH- 
Note that this spacelike singularity r = rsing is not a 
central singularity r — Q. The physical singularity, where 
K — Kpianck ) practically coincide with the mathematical 
one, where K = oo. 

Finally we note that the laws (|35() . (|36|l has been con- 
firmed by numerical calculations. 



B. Massless scalar field 

The case of a scalar field has been analysed by Burko 
in Here we extend his analysis. 

This case differs drastically from the case of dust. In 
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FIG. 2: Critical value as a function of eo for a) Dust case, b) Scalar case and c) Ultrarelativistic gas. 



this case a mathematical singularity Vging can exist only 
at Tsing ~ 0. In the vicinity of this singularity the solu- 
tion ((JSJ, (Pn|l can be written in the first approxi- 
mation as follows: 



9tt 

Qrr 

$ 



-(/? 



^ 2m 



(3+2 



VP+Tlnr, 



(37a) 
(37b) 
(37c) 



where m, C and (3 are constants. Wc note that we use 
constants m and C which are different from Burko's con- 
stants. Our constants have direct physical meanings: m 
is the black hole mass, C is a gauge parameter, related 
to the possibility of changing the scale of measurement 
of the t space coordinate. Also we have (P = ||^^^| 4m^C 
where d is a constant used by Burko in ^ . 

The exponent (3 depends on the amplitude of the scalar 
field. As Burko demonstrated j3 > Q \i q ^ Q. So in 
the vicinity of the singularity the value of ( ^) , which 
is the only term in the equations H25|) . H2t)|l . (^3011 . which 
determines the strength of the scalar field, can not be 
smaller then (unless it is equal to zero identically). 
To understand the behaviour of the singularity in this 
case let us note the following; 

In the metric (I21|l the scalar field can be represented as 
a sum of two equal fluxes moving in opposite directions 
along the (spacelike) t-axis with the fundamental velocity 
c. Indeed, let us suppose that in u, v coordinates there 
are everywhere and always two equal, opposite directed, 
fluxes along these coordinates, hence ^ — ^ which 
depend on r = u -\- v, but not ow. t = u — v. Then there 
is a coordinate transformation: 

R-t,v ^ R + t (38) 

which corresponds to a transformation to the metric 124|l 



but with another time coordinate R : dR 



The transformation H38I) corresponds to a transformation 
of the tensor of the scalar field: 



^ tt — UU I vv 



AU other Tiu = 



(39a) 
(39b) 
(39c) 
(39d) 



Applying transformation H38() to H30() . the new energy- 
momentum tensor depends on the timelike coordinate 
r but not on t. The existence of two opposite fluxes 
near the Cauchy Horizon should lead to two nonlinear 
effects: mass inflation and shrinking of the CH down 
to r = 0. The uniformity and equality of the two fluxes 
lead to the situation where both effects manifest themself 
simultaneously and there are not any gradients in space. 
To see these effects we perform the numerical integration 
of the system (I25II . (|26|I . H3()() . Here, as for the case of 
dust, we start the computation from rinu — 0.95 • teh ~ 
1.25, put the initial values of gu and equal to their 
values for the zero matter content Reissner-Nordstrom 
solution (with q = 0.95, m = 1.0) at Tinu and vary the 
characteristic of the initial amplitude of the scalar field, 
eo- 

In fig. 121 one can see the propagation (r vs. t) of the 
ingoing signal with the velocity c (c = 1) in models with 
different eo- Fig- El shows the mass function as a function 
of r for the same choices of eo and fig. [Slpresents examples 
of the evolution of the metric functions, gu and and 
K in the models with different eo- Also we refer to fig. 
121 (line b) depicting function of eo- 

From these figures it is clearly seen that in models with 
very small eo (e.g. eo = 0.0001) there is a manifestation 
of a mass inflation at r close to the CH. First of all we see 
that the light signal propaga tes al ong r w rcH during a 
long period (line 



in fig. 3(a) I 



-dr. 



This is a nessesary 
condition for mass inflation to occur. Secondly, we see 
more directly that the massfunction, which was small at 
large r, starts to manifest mass inflation at r close to 
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to = 0.005 



(b) Line a) eo = 0.010, b) eo = 0.025, c) eo = 0.050, d) 
eo = 0.250 



FIG. 3: r versus t for the scalar case for various eq. 
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(a) Line a) eo = 0.0001, b) eo = 0.001, c) eo = 0.0025 and line 
d) eo = 0.005 



(b) Line a) eo = 0.010, b) eo = 0.025, c) eo = 0.050 and line d) 
eo = 0.250 



FIG. 4: Mass function versus r for the scalar case for various eo- 



rcH (line "a" in fig. 4(a) I. The metric functions gu and 
grr behaves like the case of the pure Reissner- Nordstrom 
solution at larg er r, b ut in the vicinity of rcH they start 
to collapse (fig. 5(a) I. We also see that K de monstr ates 
a sudden sharp increases at r close to rcH (fig- |5(d)| line 



"a" ) , and it reaches the Kpianck value at Tc close to rcH 
befor e the shrinkage of the CH manifests itself strongly 
(fig. |5(d)| )- Thus here we have the physical singularity 
at r close to rcH 7^ 0. 

At larger eo, the term associated vifith scalar matter 
in the Einstein equations starts to dominate over a term 
which represents the electric charge much earlier, hence 



the manifestation of the electric field (which is respon- 
sible for the origin of the CH) is not so essential in this 
case. Functions gu and g^ differ from the case of the 
Reissner-Nordstrom at r essentially larger than rcH (see 
figs. 5(b)|5(cj and compare with fig 5(a) which essen- 
tially behaves like the Reissner-Nordstrom solution for 
r > rcH )- For the cases eo > 0-001, the light signal 
propagates at r close to rcH for a very short period of 
time (hne "b,c,d" on Fig. |3(a)| l- For eo > 0-01 the light 
signa l does not feel the presence of rcH at all (see fig. 
3(bl. For 0.01 < eo < 0.03, we observe the shrinkage of 
the model to r close to r = before K reaches Kpianck 
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(d) Kretschmann scalar for a) eo = 0.0001, b) eo = 0.01 and c) 
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FIG. 5: Metric functions (a-c) and Kretschmann scalar (d) versus r for the scalar case for various eo- 

(see line "b" on fig. So for these values of eo the C. Ultrarelativistic gas, -P = f 

physical singularity is at r close to r = 0. 



At big values of eo (for example eo > 0.050) there is not 
any manifestation of the effects near r — tch because in 
this case the light signal does not propagate long enough 
along r « tch for mass inflation to occur. The mass 
function nevertheless increases impetously with decreas- 
ing r due to the compression of the model and K reaches 
the critical value Kpianck at rather big r (see figs. |21|5(c)| 



and 5(d) I 



The case P = ^ is in some sense intermediate between 
the cases of pressureless dust and scalar field as it is seen 
in fig. 121 In fig. El one can see the contraction of the 
model and corresponding increase of the Kretschmann 
scalar K. There is not any manifestation of the mass 
infiation near r « rcH, but only the nonlinear effect of 
the increase of the mass function because of the matter 
squeeze. 
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FIG. 6: Metric functions (plot a) and Kretschmann scalar (plot b) versus r for the case of ultrarelativistic gas for eo = 0.02. 



V. PHYSICS OF THE INTERIOR 

In this section we analyse the fully nonlinear processes 
inside a spherical, charged black hole with a scalar field, 
as described in section ^ using results from numerical 
simulations. Our numerical code is described in appendix 
IXIand tested in appendix IbI As mentioned in the intro- 
duction, some parts of this problem have been discussed 
in works [Tol [13 . ITHI |. In this section we extend these 
analyses and reveal new aspects of the problem. In sub- 
section IV Al we investigate a simple compact pulse. In 
subsection IV Bl we investigate a somewhat more compli- 
cated compact pulse and in subsection lV Cl we investigate 
the influence of the T^u flux on the singularities. 

To perform this analysis we specify different boundary 
conditions along some initial u — uq = 0.00 and w = uo = 
5.00 to imitate some physical fluxes of energy into the 
charged black hole, perform numerical simulations and 
analyse the results. Throughout this section, the black 
hole, prior to influence from scalar pulses, has initial mass 
mo = 1.00 and charge q = 0.95. Also, our domain of 
integration is from 5.0 < v < 20.0 and 0.0 < u < 30.0. 
For all simulations the gauge is chosen as described in 
section lUl 



A. Simple compact pulse 

We start from the simplest case when the flux of the 
scalar field into the charged black hole is specified along 
initial u ^ uq outside of the black hole in the following 
way: 



$_„(moj — v4sin^ 



Vl - Vq 



(40) 



where Vq and Vi marks the beginning and end of the ingo- 
ing scalar pulse, respectively (i.e. we set the beginning of 
the pulse equal to the beginning of our computational do- 
main) and A measures the amplitude of the pulse. This 
can readily be integrated to give: 



A 
47r 



2tT {v - Vo) - {Vi - Vq) ! 



2tt- 



Vl 



(41) 

After the pulse, at v > vi, the flux through u = uo is set 
equal to zero, i.e. ^,v{uo,v) = 0. The flux of the scalar 
field through initial ingoing segment t; = vg is set equal 
to zero: $_„(«, ijq) = 0. This means that there is no flux 
of energy from the surface of a collapsing charged star 
into the computational domain. 

Note that we formulate the initial condition directly 
for the flux Ty^ = (<I'_„)^/47r of the scalar field through 
the surface u = uq, rather than for <I> itself since Ty^^ has 
the direct physical meaning. Also remember from section 
im that once the flux through the two initial surfaces has 
been chosen, all other initial conditions are determined 
by our choice of gauge and the constraint equations. 

In our computations we vary the width of the signal 
A — {vi — Vq), and its amplitude A, in a broad range. 
In fig. [7|is seen a typical example of the evolution of the 
scala r fie ld <I> for the case of A = 1.00, A = 0.05. Fig. 
7(a) and 7(b) represents the evolution of t he flu x Ty^oi 
the scalar energy into the black hole. Fig. |7(c)| and 
shows the T^u flux which arises as a result of Ty^ being 
scattered by the spacetime curvature. In fig 7(a)|7(b) the 
initial pulse (between 5.0 < v < 6.0) and subsequent tails 
with resonances are clearly seen. In different regions, r„ti 
and Tyy are converted into one another due to curvature 
and resonances. In some regions, Tyy is locally greater 
than Tuu, it is especially noted that the highest loc al flux 
is Tyy inside the pulse (between 5.0 < w < 6.0, flg. 7(b) I. 
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(c) Lines of constant logj^Q (Tu„). Lines are from (d) logj^Q (Tuu) along v = 10.00. 

logio {Tun) = -10.0 to logio (Tuu) = -2.00 in 
A logjQ (Tuu) = 0.25 intervals. Thick dotted line marks 
logj^Q (Tuu) = —3.25. Fully drawn thick line marks apparent 
horizon. 



FIG. 7: Tvv and Tuu for the simple compact pulse, case: A = 1.0, A — 0.05. 



We will now consider some direct effects related to 
these fluxes. 



1. Focusing effects 

One of the first noticeable effect is that the initial pulse 
leads to an initial change of the outer apparent hori- 
zon (OAH) and inner apparent horizon (lAH) in the re- 
gion within the pulse itself, e.g. 5 < u < 6 for fig [3 In 



fig. |8(a)| it can be seen that the mass function, m (eq. 
(|22|l l near the lAH increases correspondingly as the lAH 
moves from u ~ 27.12 to u « 27.51. A similar change 
can be seen for the OAH in the same region (e.g. fig. 



9(a) I. The increase of the mass function and the change 
of the apparent horizons in this region is the trivial effect 
of mass being pumped into the black hole by the Tu„-fiux 
of the initial pulse. 

The dramatic change of the lAH at ?; « 7 — 8, however, 
is related with other nonlinear effects. We remember that 
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(a) lAH and mass function. Line a) is the lAH (left and 
bottom axis). Lines b)-d) represents the mass function along 
u = 27.338, u = 27.533 and u = 27.884 respectively (right and 
bottom ajcis). 



(b) Mass function along lines of constant u. Separation 
between lines is Am = 0.40, bottom line is along u = 24.00, top 
line is along u = 30.00, thick dotted line is along u = 26.00. 



FIG. 8: Illustrations of the mass function for the simple compact pulse, case: A — 1.0, A = 0.05. 




(a) Lines are from r = 0.680 (bottom line) to r = 0.660 (top (b) Lines are from r = 0.68710 to r = 0.68665 in 

left line) in Ar = 0.001 increments. Thick dotted line is Ar = 5 ■ 10~^ increments. Thick dotted line is r = 0.68690. 

r = 0.669. 



FIG. 9: Lines of constant r for the simple compact pulse, cases (a): A = 1.0, A = 0.05 and (b): A — 1.0, A = 0.01. Fully 
drawn thick line marks apparent horizon. 



worldlines of imaginary test photons along u — const and 
V — const are under action of the gravity of the radiation 
Tyy and T„„, which leads to a focusing effect. For exam- 
ple, in the absence of scalar radiation, outgoing photons 
along u = const slightly above of the lAH will go to 
greater r as t; — s- oo in the Reissner-Nordstrom solution. 
With the existence of scalar radiation, a similar outgoing 
ray initially slightly above the lAH will now, because of 
the focusing effect of the T^y and T„„ radiation, go to 



smaller r and generate a maxima = 0) , which 

o \ dv / 7'— const 

cor respon d to the position of th e lAH . This is seen in 
fig. 9(a) and more clearly in fig. |9(b)| This effect leads 
to a drastic change of the shape of the lAH. It is seen 
by lines c) and d) in figure 8(a) that this change of the 

is not accompanied by 
About the 



lAH in this region (w « 7 - 
any significant change of the mass function 
increase of the mass function at v > 8 see below. 

In the case of smaller initial amplitude of the pulse A, 
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(a) Lines of constant logj^Q (Tuu)- Lines are from 
logio C^"") = -10.0 to logio (Tuu) = 0.00 in 
A logj^Q (Tuu) = 0.50 intervals. Thick dotted line marks 
logjg (Tuu) = —1.50 (decreasing "outwards"). Fully drawn 
thick line marks apparent horizon. 



(b) logio C^"") along v = 15.00 
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(c) Lines of constant r. Lines are from r = 0.520 (bottom loft 
line) to r = 0.475 (top right line) in intervals of Ar = 0.005. 



(d) Lines of constant logjQ (Tgg). Lines are from 
logio {Tge) = -2.0 (bottom right line) to logio {Tgg) = 32.0 
(top right line) in intervals of A logio (^se) ~ 2.0. 



FIG. 10: Various contour plots for the simple compact pulse, case: A = 1.0, A — 0.20 



the change of the shape of lAH due to focusing starts 
later. For example, for the case A — 1, A = 0.01 (see 
figure 9(b) I, the change starts at v ^ 10. In this case the 
change of the OAH and lAH in the region 5 < v < 6, 
related to the initial pulse of the scalar field, is so small 
that it is invisible in the figure. 



In figs. |10(a)|lG(c) (case: A = 1.0, A = 0.2) one can 

see another manifestation of the focusing effect related 
with the change of the flux T„„. The figures shows close 
correspondence between T„„ and the rate of focusing of 
lines of constant r. Especially we note that along the line 



u ~ 22.5 for V > 8 there is a minimum of T^ u (seen by the 
collection of very closely spaced lines in fig. 10(a) and as 
the local minima in fig. 10(b) i. Comparing with figure 
10(c) we see that the lines of r = const shows minimal 



focusing along this minima, compared to the focusing 
at u < 22.5 and 24 < u < 28. At u > 28 there is a 
decrease of Tuu and we see a corresponding decrease of 
the focusing effect. 

We should also remember that in the dynamic equa- 
tions (equations H10|l - H12|) ') and the expression for Tgg, 
(eqs. ^ + O) the scalar field appears only in the form 
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(a) Mass function along lines of constant u. From u = 22.40 
(bottom line) to u = 30.00 (top line) in Au = 0.40 intervals. 
Thick dotted line is along u = 26.00. 



(b) Lines of constant logj^Q (Tw)- Lines are from 
logio {Tw) = -10.0 to logio (T„„) = -0.25 in 
Alogio {Tw) = 0.25 intervals. Thick dotted line marks 
logj^Q {Tyv) = —4.75. Fully drawn thick line marks apparent 
horizon. Closely spaced lines near v ^ 15 inside apparent 
horizont marks local minima. 



FIG. 11: Plots for the simple compact pulse, case: A = 1.0, A — 0.10. 
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(a) Mass function along v = 20.00. 



(b) logio {Tuu) along v = 20.00. 



FIG. 12: Plots for the simple compact pulse, case: A = 2.0, A — 0.20. 



of the product ^.u^^v and the same nonhnear effects can 
be described in terms of Tgg component which is pre- 
sented in fig. 10(d) For example, where Tuu has its 
minima, a corresponding effect is clearly visible in fig. 
10(d)| We also note from fig. ^| that for this case the 
initial pulse is so strong that its flux changes the lAH 
inside the pulse and there are no double turns as it was 
seen on fig. \7\ 



2. Mass function 



Let us come now to the discussion of the behaviour 
of the mass function m. As we remember, the (modest) 
increase of m 



8(a) 



IS 



in the region 5 < w < 6 in fig 
related to the input of the energy in the initial pulse 
The increase of m seen in the region 8<y<15infig.|Sl 
is related partly with the compression-effect (see section 
IlH) . but still it is very difficult to separate this effect from 
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FIG. 13: Mass function and Kretschmann scalar along lines of constant u for the simple compact pulse, case: A = 2.0, A = 0.2. 
Separation between lines is Au — 0.40. Lines are from u — 19.40 (lowest right line) to u — 29.80 (near vertical line in lower 
left corner). Thick dotted line is along u = 24.60. 



the beginning of the mass infl ation. The essentially faster 
increase of m at w > 15 (flg. 8(b) I is the manifestiation 



of the mass inflation when we come to the CH. 

Mass inflation depends mainly on the Tyy flux but also 
on the flux and, as described in section UTTl the co- 
existence of both fluxes is essential for mass inflation to 
occur. An example of the dependence on Ty^ is seen in 
flg. ^2 ( case A = 1.0, A = 0.10). It is seen that where 
Tyy has a minimum (the narrowly spaced lines at w « 15) 
the increase of m almost stops. Fig. E| demonstrates 
the dependence of mass inflation on T^u for the stronger 
pulse: A = 2.0, A = 0.20. This pulse is so strong that a 
r = singularity is formed in the domain (this is further 
discussed in the next subsection). At the minimum of 
Tyu at u w 23.8 the increase of m also stops and at 
M > 24 where Tyu increases rapidly, m also has similar 
rapid increase. However, the line plotted terminates at 
the r = singularity, thus the flnal rapid increase of mass 
is a combination of the effects of compression and mass 
inflation. 

The compression effect can be more clearly seen in flg. 
which shows m, along lines of constant u, again 
2.0, A = 0.2. When one comes to the 



13(a) 

for the case A 
r — singularity, compression tends to inflnity and we 
see catastrophic inflnite increase of m. This can be seen 
by the near vertical lines in the lower left in the flgure, 
which represents lines of high u. These lines experience 
a catastrophic inflnite increase of mass as they approach 
r = 0, as indicated by these lines being near vertical. 

The remaining right hand side lines which run in all 
the range 5 < v < 20, represents lines of constant u 
which reach the CH and the mass increase along those 
lines are due to the mass inflation. The line marked by 



thick dashes represent a line that comes close to the point 
where the r = and CH singularities meet. The struc- 
ture of this line is more complicated as it is influenced by 
both processes. 

Finally, in flg. \iWh) 
for the same lines. 



is seen the Kretschmann scalar 



3. The singularity 

Now we will discuss the singularity. When the initial 
pulse is rather weak we cannot see the manifestation of 
the spacelike singularity in our computational domain. 
However, we can see the asymptotic approach of u = 
constant test photons to the CH singularity. In flg. 14(a) 
we see that for the case A — 2.0, A — 0.01, all our test 
photons come asymptotically to the same value at r « 
0.69, corresponding to the analytical value for r at the 
CH for the pure Reissner-Nordstrom solution, i.e. the CH 
singularity itself does not show any tendency to shrink 
down (within our computational domain). 

With an increase of the amount of energy in the initial 
pulse one can observe a nonlinear effect of shrinkage of 
the CH-singularity under the action of the gravity of the 
irradiating flux T„„ together with the Tyy flux. In flg. 
14(b)| (case A = 2.0, A = 0.125) the test photons with 
greater u — constant comes asymptotically to smaller 

(case A = 2.0, A = 0.180) 



14 c 



values of r and in flg. 
the pulse is so strong that the CH almost (but not quite) 
shrinks down to r = 



In flg. 

both the 



case A = 2.0, A = 0.200) one can see 
manifestation of the shrinkage of the CH sin- 



gularity (photons with higher u come asymptotically to 
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Case: A = 2, A = 0.010. From u = 23.00 (rightmost line) 
to u = 29.80 (leftmost line). 





to n = 29.80 (leftmost line). 




(c) Case: A = 2, A = 0.180. From u = 20.20 (rightmost Une) 
to u = 29.80 (leftmost line). 



(d) Case: A = 2, A = 0.200. From u = 19.80 (rightmost line) 
to u = 29.80 (bottom leftmost line). 



FIG. 14: V versus r along lines of constant u for amplitudes for the simple compact pulse. Separation between lines of constant 
u is Au = 0.40. 



smaller r) and existence of the r — singularity (photons 
with the highest u come to r = 0). 

Figure [T5I shows lines of constant r and the position of 
K = Kpianck (marked by the thick dotted line) for the 
two strongest cases from fig. El We remember that this 
line and places with higher K should be considered as a 
singularity from the point of view of classical physics. 
Thus, for both these cases the physical singularity is 
placed at finite values of v and is not a null singularity. In 
fig. 15(b) we furthermore see the r = spacelike singu- 
larity inside of our computational domain. This singular- 
ity can be considered as a result of mutual gravitational 
focusing of Tyy and T^u fluxes in the region between in- 
ner and outer apparent horizons. At small w < 15 the 



physical singularity practically coincide with r = 0, but 
for i; > 17 we see that the spacetime structure of the 
physical singularity is quite different from the structure 
of the mathematical singularity. The physical singularity 
here depends mainly on the true CH-singularity, but its 
position in the u — v diagram is quite different from the 
position of the true CH-singularity which is at w = 00. 



This can also be compared with fig. 13(b) from which 
we see different behaviours of K for the test photons for 
the strong case of A = 2.0, A = 0.200. The fines in the 
lower left hand corner, sharply increasing to near vertical, 
are the lines which come to r = 0. The thick dotted line 
is the line which comes to a point at the singularity close 
to the meeting of the r = and CH singularities. The 
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(a) Case:A = 2.0, A = 0.18. Lines are from r = 0.45 (bottom (b) Case:A = 2.0, A = 0.20. Lines are from r = 0.42 (bottom 

left line) to r = 0.03 (top right line) in intervals of Ar = 0.01 left line) to r = 0.20 (top right line) in intervals of Ar = 0.01. 

Top right thick line mark r = singularity. 



FIG. 15: Lines of constant r for two different simple compact pulses. Thick dotted lines in right part of the figures represents 
line of planckian curvature. Bottom fully drawn thick line marks OAH. 




(a) Lines of constant logj^Q (T„i,). Lines are from (b) Lines of constant r. From r = 3.00 (bottom line) to 

logj^Q {Tvv) = —10.0 to logj^Q (Tvv) = 2.0 in intervals of r = 0.10 (top right thin line) in Ar = 0.10 intervals. Thick 

^l°gio C^vv) = 0.50. dotted line marks K = Kpia„ck- 



FIG. 16: Contour lines for double sine pulse of form of eq. (I42I I. case A = 2.0, A = 0.25. Thick bottom line is AH, thick upper 
line is r = 0. 



remaining lines in the right hand side are the lines which B. Double sine pulse 

come to the CH-singularity. 

In this subsection we choose the ingoing flux to be of 
the form: 



^^y{uo,v) = VgA cos ( TT 







\ ■ ( w - «o \ 






1 sm TT 









(42) 
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(a) Lines of constant logj^Q (Tuu)- Lines are from 
logio (Tun) = -10.0 to logio (Tuu) = 3.00 in 
Alog^o (Tuu) = 0.25 intervals. Thick dotted line marks 
logj^Q {Tuu) = 1.50, decreasing leftwards. 



(b) Lines of constant logj^Q (Tw)- Lines are from 
logio C^vv) = -10.0 to logio (Tvv) = 3.00 in 
Alogj^Q (Tuv) = 0.25 intervals. Thick dotted line marks 
logio {T^vv) = 1.50, decreasing downwards. 



FIG. 17; Contours for the double-flux case based on the simple compact pulse: A = 2.0, A — 0.25. Thick top right line marks 
r = 0. 
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FIG. 18: Lines of constant r for the double-flux case based 
on the simple compact pulse: A — 2.0, A — 0.25. Thick top 
right line marks r = 0. Lines are from r — 2.50 (bottom left 
line) to r = 0.10 (top right thin line) in intervals of Ar = 0.10. 
Thick dotted line marks K = Kpianck- 

instead of equation H4(J|I . This readily integrates to give 
the initial expression for ^(uq,v): 




where vq and vi as before, marks the beginning and end 
of the pulse respectively and A is the amplitude of the 
pulse. Also, as before we set ^^v{uq,v) — for v > vi. 
The pulse is scaled in such a way that for a given width 



A — vi ^ vq and amplitude A, the integral of the initial 
flux, J_^^ Tyy dv, is equal for pulses of the form l|40() and 
(|42ll . Eq. (|42|l has the shape of a double pulse, rather 
then H4Uf) which is the shape of a single initial pulse. 
This pulse is more complicated than (|40|l . but is similar 
in shape to the pulse shapes used in some other papers 

(e.g. mM)- 

We have performed a series of computations based on 
eq. (|42|l (all other initial conditions equal to those in the 
previous subsection). The results of these computations 
demonstrate a more complex picture of interplay between 
the scalar fluxes than in the case of subsection lV Al This 
is natural because of the more complex shape of the ini- 
tial pulse. But the main physics and principal properties 
of the singularities are the same in the two cases. An ex- 
ample illustrating the increased complexity in structure 
of the fluxes can be seen in figure^] which illustrates the 
flux for the case:A = 2.0, A = 0.25. The shapes of 
the apparent horizons and the central singularity r = 
are now more complex as well as the distribution of the 
Tyy field. Still the general characters of the r = and 
the physical K ~ Kpianck singularities are the same. 



C. The influence of the Tuu flux 

So far in all our analyses we have assumed that the 
flux of the scalar field T„„ through v — vq was zero and 
that any Tun-flux arised only as a result of scattering of 
the T^^-flux by the curvature of the spacetime. Now we 
would like to consider the influence of a T„„-flux through 
the surface v = vq. To do this we will consider an ex- 
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(a) Mass function along lines of constant v. Lines are from 
= 7.04 (bottom right) to v = 19.04 (top left) in Ad = 1.0 
intervals. 



(b) Mass function along lines of constant u, corre spondi ng to 
the symmetrical equivalent of the lines in fig. |19(a)| 



FIG. 19: Mass function for the double-flux case double-flux case based on the simple compact pulse; A = 2.0, A = 0.25 
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(a) Kretschmann scalar along lines of constant v. Lines are 

from V = 7.04 (bottom right) to v = 19.04 (top left) in 
Ati = 1.0 intervals. Horizontal dashed line is K = Kpi^nck- 



(b) Kretschmann scalar along lines of constant u, 
corresponding to the symm etrical equivalent of the lines in fig. 



FIG. 20: Kretschmann scalar for the double-flux case double-flux case based on the simple compact pulse: A = 2.0, A = 0.25 



treme case when T^u through v = is equal exactly to 
T^jy through u — const just inside of a black hole. More 
concretely we do the following; We specify some initial 
^(uo,v) along u = uq with a pulse width A = vi — vq 
and amplitude A and set all initial conditions equal to 
those in subsection IV Al including $.„(it, uq) = along 
V = vo- These initial data are then simulated as usual, 
however only in the computational domain of 5 < w < 20 
and < u < uahi where uah is the first computa- 
tional point which is inside of the outer apparent hori- 
zon along V = vi. Because of scattering of the initial 



pulse there is now a T^y flux into the black hole along 
the line u = uah for v > vi. We then stop the com- 
putation and start a new with the following domain of 
integration: vi < v < 20 and uah ^ u < Umax (where 
Umax - UAH = (20 - -wi)). Aloug the (new) outgoing 
initial hypersurface, u = uah, all the variables are kept 
as they were in the original simulation, while along the 
(new) ingoing initial hypersurface, v = vi, initial data 
are set equal to the data along the outgoing hypersurface, 
hence we have completely symmetrical initial conditions. 
Subsequently we performed a computation for the new 
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(a) u vs. r along lines of constant v. Lines are from v = 7.04 
(rightmost) to ti = 19.04 (bottom left) in Av = 1.0 intervals 



(b) V vs. r along lines of constant u, correspo nding t o the 
symmetrical equivalent of the lines in fig. |21(a)| 



FIG. 21: u and v versus r for the double-flux case double-flux case based on the simple compact pulse: A = 2.0, A = 0.25 



computational domain. It is obvious that all conditions 
in this domain are symmetrical with respect to u and v 
and hence the boundary fluxes T^u along v — vi and 
along u — UAH must be exactly equal. 

We performed our computations for the compact sim- 
ple pulses (eq. if^ ) for different parameters A and A. 
An example of our results can be seen in figs. ll7l-l21lfor 
the case A = 2.0,^ = 0.25. All pictures are symmetrical 
with respect to u and v, as they should be. The outer ap- 
parent horizon is naturally outside of our computational 
domain. Figs. 1171 and 1 181 shows the T^u and distri- 
butions and r contour lines. On these figures there is no 
inner apparent horizon because it coincide with another 
CH singularity at m ^ oo (see below). But on fig. [T5I 
one can see the mathematical strong singularity r — Q 
and the physical singularity (with K = Kpianck) which 
is reached long before the CH singularity at u ^ oo and 
?; ^ oo in most of the computational domain. 

Fig. ^1 shows the huge increase of the mass function 
with growing u and v, and fig. 1201 shows the increase of 
the Kretschmann scalar as function of u and w up to if = 
Kpianck- Finally fig. 071 shows light signals along v = 
const and u — const respectively. This figure shows both 
the existence of null singularities of the CH singularity- 
type and r = singularity. Indeed for example in the 
case fig. 21(a) one can see that for some v = const (the 
rightmost curves) for u — s- oo, signals goes asymptotically 
to practically r = const (asymptotically approach to the 
null singularity at m — > cxd). Smaller asymptotic values r 
are seen for the signals at bigger v corresponding to the 
focusing effect. Finally the signals with v big enough, 
come to the central singularity r = 0. The symmetrical 
picture for the signals with u = const is shown in fig. 
21(b) Here we see the asymptotic approach to another 
null singularity at f ^ oo. 



Note also the different character of the increase of mass 
function on fig. 1191 for lines of big and small v and u re- 
spectively. For example on fig. 19(a) the topmost lines 
(big V — const) come to the central singularity but lines 
with small v — const (bottom lines) go to the null singu- 
larity at w — > oo. The different behaviour of these lines 
correspond to the different nonlinear processes influenc- 
ing the mass function for lines terminating at the r = 
singularity and reaching the CH singularity. 

Thus in this case there are three different singularities 
inside the black hole: two null singularities at u — > cxd and 
V —^ oo and the physical singularity at K = Kpianck- In 
addition there is also a central (mathematical) singularity 
at r = 0. 



VI. CONCLUSIONS 

In this paper we investigated the physics of nonlinear 
processes inside of the spherical charged black hole, non- 
linearly perturbed by a selfgravitating, minimally cou- 
pled, massless scalar field. For this purpose we created 
and tested a numerical code which is stable and second- 
order accurate. For our computations we used an adap- 
tive mesh refinement approach in ingoing it-direction. 

The following nonlinear physical processes are impor- 
tant inside the black hole: Scattering of radiation by the 
curvature of the spacetime, gravitational focusing effect, 
mass inflation and squeeze of matter with pressure. 

At the beginning of our analysis we used a homoge- 
neous approximation to clarify some physical processes 
near a spacelike singularity. In this approximation one 
supposes that near the singularity temporal gradients are 
much higher than the spatial gradients, so one assumes 
that all processes depend on the time coordinate only 
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(uniform approximation) . We used both analytical anal- 
yses and a numerical approach to analyse three different 
physical matter contents: dust, a massless scalar field 
and an ultrarelativistic gas. 

For the case P = (dust) we found that the singularity 
is at r = Tsing ^ 0, rcH < Tsing < teh, where tch 
and veh sltc the positions of the Cauchy Horizon and 
the Event Horizon when there are no dust at all. The 
value Tsing decreases monotonically with a decrease of 
the matter contents and tends to Vsmg = rcH when the 
matter contents goes to zero. 

In the case of the scalar field, the uniform approxi- 
mation demonstrates more a complex behaviour. Here 
the scalar field can be represented as a sum of two equal 
fluxes moving in opposite directions. One can for this 
case see the manifestation of both the effect of mass in- 
flation and the effect of shrinkage of the CH down to 
r = 0. For very small matter contents (eo <C 0.01) the 
Kretschmann scalar, K, becomes equal to the Planck 
value at r close to rcH- So in this case the physical sin- 
gularity (when K = Kpianck) is at r « rcH- For larger 
values of eo, (e.g. 0.01 < eo ^ 0.03), the CH does not 
manifest itself and the model squeezes to r very close to 
r = before K reaches Kpianck- It means that for these 
values of eg, the physical singularity practically coincides 
with r = 0. For rather big eg, K reaches Kpianck at 
rather big r as it was for the case of dust, P = 0. 

In the case of matter with isotropic relativistic pres- 
sure, P = e/3, we have the situation intermediate be- 
tween P — and the scalar field. The physical singular- 
ity, in this case, is located at r essentially greater then 
r = 0. 

We performed the analysis of the full nonlinear pro- 
cesses inside the spherical charged black hole with a 
scalar field using the numerical approach. This analysis 
extends the analysis of the earlier works ,1^ ^1^ JJj and 
reveal new aspects of the problem. The detailed descrip- 
tion of the results is given in sectional We analysed non- 
linear gravitational interaction of the fluxes of the scalar 
field, the dependence of the effects on the boundary con- 
ditions, analysed the focusing effects, mass inflation and 
squeeze effect and the behaviours of the Kretschmann 
scalar K. We payed special attention to the analysis 
of the singularity in subsection IV A 31 We investigated 
the focusing of the CH singularity and its dependence on 
the boundary conditions. We determined the position of 
the physical singularity (where K = Kpianck) inside the 
black hole and demonstrated that this position is quite 
different from the positions of the mathematical r = 
singularity and CH singularity. 

The results mentioned above were obtained with the 
scalar flux into the black hole in the form of a simple 
compact sine-pulse with different amplitudes and widths. 

We also analysed the physics in the case of a scalar flux 
in the form of a double sine-pulse qualitatively similar 
to the usage in [2^ In this case physics is more 

complicated, but the main characteristics of the results 
are the same as for the simple pulse. 



Finally we investigated the influence of the boundary 
T-uu-flux on the physics of the singularity. We demon- 
strated that it is possible to have the existence of three 
different singularities inside the black hole: two null sin- 
gularities at u — > cxD and v —^ oo and the physical singu- 
larity K — Kpianck in addition to a central mathematical 
singularity r = 0. 
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APPENDIX A: THE NUMERICAL CODE 

In this appendix we describe the structure of the nu- 
merical code used to obtain the numerical results pre- 
sented in section Ivl 

1. The numerical scheme 

Our approach to numerically integra ting the field equa- 
tions is similar to that of other authors [ltt[l4llalla . l28l| . 
Denoting the three unknowns, r, a, $ as hj with j = 1, 3, 
the three evolution equations (|10|I - H12I) are each of the 
form: 

hj,uv = Fjihk,hk,v,hk,u),k = 1,3 (Al) 

We denote four grid points in the (u, ?;)-domain by pi = 
{u,v), p2 = {u,v + Av), = {u + Au,v) and p4 = 
(u + Au, v + Av), where Av and Au are finite increments 
in the v and u directions (see fig. I22|l . 

By combining a Taylor expansion hj at these four 
points (evaluated around an intermediate point p^ = 
{u + Au/2, V + Aw/2)) with eq. (IJTl we find that: 

hj(p4) = hjips) + hj{p2) - hj{pi) 

+AuAv[Fj{hk ip5 ),hk,u{P5),hk,v{Pb)) 
Au ^ Av^ 

+0{Au^,Av%k^l,3 (A2) 

i.e. we can evaluate hj at p^ with second order accu- 
racy if we know hj at the four points pi,p2,P3,P5 and 
hj,v , hj^u at P5 . Initially we only know the values of hj 
along in- and outgoing null segments (see section ITU) i.e. 
for example at pi,P2 and p^. However, by employing eq. 
IA2I twice in a predictor-corrector style we can evaluate 
hj , hj^y , hj^u at P5 with the desired accuracy and hence 
subsequently find hj at p^. 
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FIG. 22: Schematic of the points used in the numerical 
scheme. 



Equation (jA2p constitutes the heart of our code. As 
described in section m the values of the unknown func- 
tions along the initial hypersurfaces uo,vq are given by 
the constraint equations © and © . We can then use our 
numerical scheme to calculate the values of hj along the 
ingoing ray v ~ vq+Av starting at the point u = uq + Au, 
then solving for u = uq + 2Au etc. until we reach the 
last grid point at u = u final- Using the (now known) 
values of hj along u = uq + Av, we can then calculate the 
unknowns along the next ingoing ray along w = wq + 2Au 
and so on throughout our computational domain until we 
reach the end of our computational domain ed v — v final ■ 



2. Splitting algorithm 

When outside the event horizon (EH), the scheme de- 
scribed above works reasonably good as it is, but when 
the domain of integration crosses the EH, the unknown 
functions exhibits extremely great gradients. This can 
be seen by considering two outgoing rays, one just inside 
and one just outside of the EH. One ray is destined to be 
trapped inside the EH, while the other must escape to 
infinity, i.e. regardless of how close the two rays were ini- 
tially, their distance will diverge as their advanced time v 
grows. Integration across such great gradients using fixed 
increments in u and v would cause the numerical errors 
to explode unless the initial increments were chosen to be 
unrealistically small. To overcome this difficulty an adap- 
tive mesh refinement (AMR) approach is crucial for an 
accurate integration of the field equations. However, the 
strong gradients are mainly in the ingoing w-direction, 
hence AMR is most important in this direction. While 
AMR in both v and u directions would be desirable, im- 
plementation of AMR in only one direction, while much 



less complicated, still produces good results. This is sup- 
ported by other authors who also employ uni-directional 

AMR Hi m. 

To determine whether AMR is needed and a cell should 
be split (or possibly desplitted), we focus on the trun- 
cation error in the numerical scheme. As can be seen 
from equation (|A2|) , the primary error introduced by the 
scheme comes mainly from the terms: 



eT(P5) 



Au^ 

~2r 



Av' 
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^J^VUUU 



(P5) 



bs) 



(A3a) 



(A3b) 



Since we use AMR only in the ingoing u direction, 
thus keeping Aw constant, we can only control £"(^5) 
by splitting, hence we will focus on this term. By in- 
troducing four new points: = (u — Au,v),p-i = 
{u—Au,v+Av),ps = {u—2Au,v),pq = {u—2Au,v+Av), 
(see fig. 1^ . we can estimate hj^^^^u (and hence e") at 
P5 with the following finite difference operator: 



^],vuuu [Pb] = AvAu \^3hj (pi) - 3hj (ps) 

-3hj (pe) + Shj (pr) - hj (^2) + hj (pi) 
{ps) - {P9)] + O {Au) (A4) 



In practice we proceed as follows; Following a calcu- 
lation of hj at p4 we use eq. (|A3a(l + ()A4(I to estimate 
the error £"(^5) involved in the calculation and compare 
it to Fj{hk{p5), hk^uiPb), hk^viPb)) in eq. (|X2|) for each of 
the variables, j = 1,3. We now require that the relative 
error should be within some fixed interval: 



Sdesplit ^ 



3 

Fjihkipb), hk,u{Pb), hk,v{Pb)) 



< SspHt (A5) 



If the relative error is greater then the threshold param- 
eter S split we perform a split. Conversely, if the relative 
error is less than the threshold parameter SdespUt a de- 
split is made. Note that because we need knowledge of 
the points pq to pg to calculate the relative error, we 
don't use splitting for the first two cells u < Uq + 2Au on 
a given ingoing ray. 

If a split is required a new point is introduced at p'^ = 
{u -t- Au/2, v). The values of hj at pg are then found by 
a four-point interpolation scheme (using data from the 
two points above and the two below pg). The numerical 
scheme is then used to calculate hj at P4 = (u -I- Au/2, v + 
Av) using the points pi,P2 and pg. If a desplitting is 
needed the reverse process takes place, i.e. pa is deleted 
and the point above is used in its place. 

To obtain the results presented in section the split- 
ting thresholds were set to Sspin = 10^^, Sdespiu = 10~^. 
The resolutions along the initial outgoing null segment 
were set to Au = 1/100 and Aw = 1/12800. The initial 
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(a) Line a) VivC?"), b) iPn{^) and line c) ipN{cr). Line d) 
shows a line with a slope of minus two. 
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(b) Line a) xCC'i) ^■iid b) x(C2)- Line c) shows a line with a 
slope of minus two 



FIG. 23: Demonstration of the convergence (without AMR) of r, 3>, cr, C\ and C2 along an outgoing ray u = 23.00. 
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FIG. 24: Ci (fig. a) and C2 (fig- b) along an outgoing ray u = 23.00 for resolutions; line a) iV = 100, h) N = 200, c) TV = 400 
and line d) = 800. 



value of Am can be set rather low, since the splitting al- 
gorithm will assure that Au always have an appropriate 
value. The value of Av on the other hand, is constant 
throughout the computational domain and must initially 
be chosen to have a suitable value everywhere. The high 
value of Av has been chosen such that Av is as compara- 
ble to Am in large parts of the interior of the black hole. 
This has been verified by numerical tests to give the best 
results. 



APPENDIX B: ANALYSIS OF THE CODE 

We have tested the code and found it to be stable and 
second-order accurate. Here we presents tests to demon- 
strate this. 



1. Basic convergence tests 



We test the basic convergence properties of the code 
by performing simulations with the same initial condi- 
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tions, but with varying resolutions and with the spht- 
ting algorithm disabled. As initial conditions we use a 
Reissner-Nordstrom spacetime (with initial mass mo = 1 
and charge q = 0.95) which is perturbed by a massless 
scalar field. We set the specific form of the scalar field 
along the initial outgoing null segment in the interval 
I'o ^ ''^ ^ ^^1 to be: 



1.5 



A_ 

in 



27r {v ~ vq) — (vi — vo) sin 2tt 



Vo 



Vl 



(Bl) 

with A = 0.05 and vq = 5, wi = 6. For v > vi we set 
$ constant to $(uo,w) — ^{uo,vi). Along the ingoing 
null segment v — vq we specify ^^u{u,vo) = (these 
initial conditions are the same as those used in section 
IVj) . As domain of integration we used 5 < t; < 100 and 
< u < 23.05 and with the gauge choice described in 
section m 

We make simulations with varying resolutions N = 
1/Au = 1/Av and examine data along an outgoing ray 
u — 23.00 which is just outside of the event horizon. 
Along this ray we calculate the average absolute differ- 
ence at each data point between the unknown functions 
for two resolutions along the outgoing ray. We define: 
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El 



<^2N\ 



(B2) 



where x)^ denotes the value of variable x at the I'th data 
point (out of a total of n points) from a simulation with 
resolution N. I.e. we basically calculate the normalized 
Li-norm of the difference between a simulation with res- 
olution N and one with resolution 2N at each point for 
the three unknowns, which is similar to the convergence 
tests used in i2a. 



In figure 23(a) are shown 7/'Ar(r), ■i/;jv($) and ■0Ar(f) for 
N = 50, 100, 200, 400, 800. We see that all the lines have 
a slope of approximately minus two (compare with line 
'd'), indicating that the functions converge with second 
order. However, we also need to show that the converg- 
ing solution is a physical solution of our field equations. 
We have tested the code against the Reissner-Nordstrom 
and Schwarzschild solutions and reproduced all known 
features (e.g. location of the event horizon as a function 
of initial mass etc.) of these spacetimes. If a scalar field 
is present there are no suitable analytic solutions against 
which we can compare the numerical results, but we can 
use the constraint equations lO and lO to test that our 
code is converging to a physical solution. 

Denoting the left hand side of the constraint equations 
© and © by C\ and Ci respectively, we calculate the 
average absolute value of C\ and C2 at each point: 
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defined simila r to eq . (jB2p . with x being either C\ or 
Ci- In figure 23(b) are shown x(Ci) and x(C'2) along 
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FIG. 25: Mass function as a function of v for the resolutions; 
line &) N = 200, b) iV = 400, c) TV = 800 and d) TV = 1600. 



u = 23.00. It is noted that x(C2) is much greater then 
x(Ci). This is mainly because the line u = 23.00 is 
very close to the event horizon and the gradients along 
ingoing rays are here very large. This results in great 
inaccuracy for coarse resolutions when calculating the 
derivatives used to calculate Ci. Related to this is of 
course the fact that we are not using AMR for these ba- 
sic convergence tests and hence the large average value 
for Ci is also an indication that splitting would be most 
desirable in this region. However, x(C'2) is decreasing 
with increasing resolution, which is the important result 
in this context. We note also that the typical absolute 
value of the biggest term in eq. lO is approximately ten 
times greater than C2 even for rough resolutions. This 
means that the constraint equation is roughly satisfied for 
this case even though C2 has very large values. For the 
highest resolution (TV = 800), x{C\) shows indications of 
a slightly decreasing convergence rate. This can also be 
traced to roundoff errors when calculating the ri,„-term 
ofCi. 

Figure 24(a)| shows C\(v) for the resolutions TV = 
100,200,400,800. It is seen that Ci is clearly converg- 
ing to zero for increasing resolution. It is also clear that 
numerical noise play a dominant role for TV = 800, espe- 
cially at large w, causing the me ntioned decreasing con- 
vergence rate effect in fig. |23(b) As noted, this numeri- 
cal noise can be traced to the numerical calculation of the 
Ttjtj-term in eq. The reason that this particular term 
exhibits large roundoff errors is mainly that r^^, k, along 
the line at which we look. This results in the subtrac- 



tion of two nearly equal numbers which results in large 
roundoff errors. Figure 24(b) shows the convergence for 
C2 again for TV = 100,200,400,800. Clearly Ci is also 
converging to a zero value, but for a given resolution, Ci 
is much larger then C\ due to the proximity to the event 
horizon and the absence of AMR in these simulations as 
explained above. 

It is also noted that although C\ and Ci seem to be 



25 



0.01 



X 0.001 



0.0001 



1e-05 



1e-06 




'split 




1e-08 



''split 



(a) a: x(Ci) and b: x(C2). 



(b) a: «/'(r), b: and c: ^{rr) 



FIG. 26: Plots illustrating convergence (with AMR) as function of splitting criterias along outgoing ray at m = 24.60. 
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FIG. 27: "I> and a as functions of r along outgoing ray at u = 24.60. Error tolerance intervals are for line a) Ssput = 10 ^, b) 
Sspiit = 10~^, c) S split = 10"'' and d) Ssput = 10"'^. In all cases Sdesput = 10~^ x Saput- 



converging with second order accuracy, the convergence 
rates of Ci and C2 are not by themself as important as 
the convergence rates of the unknown functions, since Ci 
and C2 are functions of the derivatives of the unknown 
function. The unknown functions, however, directly mea- 
sures the convergence rate of the code as these are the 
actual variables being evolved. 



2. Convergence with AMR 

Another test to show that the numerical solution is 
converging to a physical solution is to monitor whether 
the code is mass conserving. Figure shows the mass 
function \ib\ along the outgoing ray u — 23.00 in 
the interval 5 < w < 100 for the resolutions N = 
200, 400, 800, 1600. We see that the mass function is con- 
verging to a constant value for high resolutions, further 
indicating that the code is converging to a physical solu- 
tion of the field equations. 

In the previous subsection we verified that with the 
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fixed increments Au = Av, our code converges to a phys- 
ical solution with second order accuracy. In this subsec- 
tion we will demonstrate the convergence properties of 
our code with the splitting algorithm turned on, as a 
function of the splitting critcrias used. 

We wish to confirm that our splitting algorithm works 
as desired and that the numerical solution converges to a 
physical solution. For this purpose we make a number of 
simulations with identical initial conditions, but vary the 
splitting thresholds. We use the same pulse shape as in 
the previous subsection (eq. jlBljl ). but with a wider and 
stronger pulse with A = 0.20 and vo = 5.0 and vi — 7.0 
which is strong enough to generate a r = singularity 
within out computational domain. Our computational 
domain is in the range 5 < w < 20 and < m < 25. We 
examine an outgoing ray {u — 24.60) which is inside the 
event horizon and in fact comes very close to the r — 
singularity. 

To measure the convergence rates we use eq. (jB2|) and 
(|B3|I only in this subsection they are functions of Ssput 
(see section instead of the resolution A''. 

Figure [26 (a) shows x(Ci) and x(C2) as functions of 
the splitting threshold SspUt (see section lA 2|l that we 
require the relative numerical error to be below for all 
computational cells, i.e. the left hand side of the figure 
denotes a more strict splitting policy. The threshold for 
desplitting is here and throughout the paper set to two 
orders of magnitude lower then the splitting threshold, 
SdepUt = ■ Sspiit- We see that as the splitting thresh- 
old Sspiit is decreased, so is the average absolute of Ci 
and C2, indicating that the numerical solution converges 
to a physical solution with decreasing error tolerance. 



For the lowest splitting threshold at the left side of 
the figure we see that lower splitting thresholds does not 
result in significant smaller values of Ci and C2. This 
arises because Ci and C2 becomes so small that numeri- 
cal noise (which is little affected by resolution) begins to 
dominate the calculations of Ci and C2. 



Figure 26(b) illustrates the convergence of the func- 
tions r, $ and a along the same outgoing ray u = 24.60. 



The figure, which is qualitatively similar to figure 23(a) 
shows the difference between simulations (V') with differ- 
ent the splittings thresholds SspUt = 10"'^, 10~^, 10"^, 
10^**, lO^'^, 10^^, 10^^. Since i/; in this subsection is a 
function of Ssput instead of N, a point in fig. 26(b) illus- 
trates the average absolute difference between the sim- 
ulation a Sspiit and the simulation with the next higher 
splittings threshold. For example the leftmost points in 
fig. 26(b) (at Sspiit = 10~^) illustrates the difference be- 
tween simulations with Sspiit = 10^^ and Ssniit — 10^^. 
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As the splitting thresholds decreases we expect to see 
-0 convergece to zero. We see that for smaller Ssput, the V' 
decrease indication a convergence of the unknown func- 
tions. This convergence is further demonstrated in figure 
[TtI which shows $ and a as functions of r close to the 
r = singularity (in the interval 0.02 < r < 0.0201) 
along the same outgoing ray that was used for previ- 
ous plots. From these figures it is clear that even in 
the vicinity of the r — singularity, the unknown func- 
tions are converging. It is also important to notice that 
even though Ci and C2 cease to decrease at rather high 
splitting thresholds, the unknown functions themselves 
continue to benefit from lower splitting thresholds. 
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